import matplotlib.pyplot as plt
import cobra
from cobra.io import validate_sbml_model
import importlib
from xml.etree import ElementTree
import utils.Model_correction as mc
import sys
import utils.model_maj as mj
import utils.viz_utils as vu
#import cplex
import cbmpy
import pandas as pd
import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
from tqdm import tqdm
from itertools import cycle
import seaborn as sns
INFO: No xlwt module available, Excel spreadsheet creation disabled CBGLPK based on swiglpk: not all methods implimented yet! ***** Using CPLEX ***** INFO: No xlrd module available, Excel spreadsheet reading disabled *********************************************************************** * Welcome to CBMPy (0.8.2) - PySCeS Constraint Based Modelling * * http://cbmpy.sourceforge.net * * Copyright(C) Brett G. Olivier 2014 - 2020 * * Systems Biology Lab, Vrije Universiteit Amsterdam * * Amsterdam, The Netherlands * * CBMPy is developed as part of the BeBasic MetaToolKit Project * * Distributed under the GNU GPL v 3.0 licence, see * * LICENCE (supplied with this release) for details * ***********************************************************************
importlib.reload(vu)
importlib.reload(mc)
<module 'utils.Model_correction' from '/home/jerem/Documents/master/RBA_cutsets/code/utils/Model_correction.py'>
HepG2, errors = validate_sbml_model("../models_storage/HepG2_v8.xml")
iHep, errors = validate_sbml_model("../models_storage/iHep_v3_genenames.xml")
HepG2_C = HepG2.copy()
HepG2.reactions.EX_m01965x.bounds = (-1000.0,-1000.0) #Glucose exchange, set to the maximal input value observed with FVA
#iHep.reactions.EX_m01965x.bounds = (-1000.0, -1000.0)
iHep.reactions.HMR_4281.bounds = (0.0,0.0) # Blocking peroxysome fermentation in healthy model
HepG2.reactions.HMR_3883.bounds = (0.0,0.0) #Serine --> Pyruvate exchange
HepG2.reactions.EX_m02630x.bounds = (0.0,0.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
HepG2.reactions.HMR_4316.bounds = (-1000.0,1000.0) # Glucose --> D-Glucitol
HepG2.reactions.EX_m02896x.bounds = (0.0,1000.0) # Serine intake
HepG2.reactions.EX_m01682x.bounds = (-1000.0,0.0) # Glucitol secretion blocked. Kept the intake just in case.
HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
HepG2.reactions.HMR_4300.bounds = (494.6155049539209,494.6155049539209)
#HepG2.reactions.HMR_4281.bounds = (100.0,1000.0) #Fermentation perox.
#HepG2.reactions.HMR_4930.bounds = (-1000.0,1000.0) # Pyruvate transfer from cytoplasm to peroxysome
#HepG2.reactions.HMR_4281.bounds = (0.0,0.0) # Fermentation in peroxysome
iHep.solver = "cplex"
iHep.objective = "biomass_components"
sol_iHep = iHep.optimize()
HepG2.solver = "cplex"
HepG2.objective = "biomass_components"
sol_G2 = HepG2.optimize()
HepG2_C.solver = "cplex"
HepG2_C.objective = "biomass_components"
sol_G2_C= HepG2_C.optimize()
df_Ci, df_Ri, df_Si, df_Mi, df_Pi, df_Xi, df_Li, df_Gi, df_Ni = vu.build_reaction_df(iHep)
df_C, df_R, df_S, df_M, df_P, df_X, df_L, df_G, df_N = vu.build_reaction_df(HepG2)
subS_HepG2 = vu.get_subsystem_fluxes(vu.build_reaction_df(HepG2))
subS_iHep = vu.get_subsystem_fluxes(vu.build_reaction_df(iHep))
subS_G2C = vu.get_subsystem_fluxes(vu.build_reaction_df(HepG2_C))
df_both = pd.concat([subS_HepG2, subS_G2C, subS_iHep], axis = 1)
df_both.columns = ["HepG2","HepG2_copy", "iHep"]
#df_both = df_both.loc[((df_both["iHep"]/df_both["HepG2"] >= 2.0) | (df_both["HepG2"]/df_both["iHep"] >= 2.0)) & ((df_both["HepG2"]>=1000.0) | (df_both["iHep"]>=1000.0))] # --> Filtre relatif + absolu
df_both = df_both.loc[(df_both["HepG2"]>=1000.0) | (df_both["iHep"]>=1000.0)] # --> Filtre absolu
clustermap = sns.clustermap(df_both.fillna(0), figsize =(10,20))
Heatmap : Représentation sous forme d'une heatmap des flux de chaque sous-systèmes actif dans trois modèles de cellule du foie différentes :
On s'attend à observer une distribution des flux dirigée vers la production d'énergie, de nucléotides, d'acides aminés ou d'acides gras pour les modèles cancéreux.
La glycolyse est effectivement plus active dans le modèle HepG2 corrigé. \ On observe également une activité importante du métabolisme du folate dans le modèle HepG2 (métabolisme des nucléotides)\ Cependant, l'activité du métabolisme des pentoses phosphates reste bien trop faible, il en va de même pour le métabolisme du pyruvate.
On s'attend à observer les différences suivantes dans le modèle cancéreux par rapport au modèle sain :
--> Résultats :
Comparaison des flux attribués à chaque groupe de réactions entre le modèle sain et le modèle cancéreux
r = HepG2.reactions.HMR_0162
r.subsystem
'Carnitine shuttle (mitochondrial)'
barplots_dict = vu.subsystem_barplots(iHep, HepG2_C, 'iHep', 'HepG2_copy',{'glycolysis / gluconeogenesis' : (10,14), 'pyruvate metabolism' : (7,11), 'pentose phosphate pathway' : (5,7), 'folate metabolism' : (5,7), 'carnitine shuttle' : (5,7), 'beta oxidation' : (5,7), "nucleotide metabolism" : (10,14)})
barplots_dict_2 = vu.subsystem_barplots(iHep, HepG2, 'iHep', 'HepG2',{'glycolysis / gluconeogenesis' : (10,14), 'pyruvate metabolism' : (7,11), 'pentose phosphate pathway' : (5,7), 'folate metabolism' : (5,7), 'carnitine shuttle' : (5,7), 'beta oxidation' : (5,7), "nucleotide metabolism" : (10,14)})
On s'attend globalement à une activité plus importante du métabolisme des nucléotides dans le modèle cancéreux. Cependant, on n'observe une activité plus importante des protéines cancéreuses que lorsqu'elles sont impliquées dans la synthèse d'ADP.
Cela pourrait indiquer un potentiel biais : l'ADP est indiqué comme produit dans la fonction de biomasse données en objectif au modèle cancéreux.
barplots_dict_2["nucleotide metabolism"][0]
barplots_dict["nucleotide metabolism"][0]
--> On s'attend à observer :
--> On observe :
barplots_dict_2["glycolysis / gluconeogenesis"][0]
barplots_dict['glycolysis / gluconeogenesis'][0]
Cette figure sert surtout en conjonction avec la précédente où on peut voir que les flux associés aux gènes PKLM et PKR sont les mêmes entre le modèle sain et cancéreux.
barplots_dict_2["pyruvate metabolism"][0]
barplots_dict["pyruvate metabolism"][0]
On s'attend à une activité de la voie des pentoses phosphates plus importante chez la cellule cancéreuse que chez la cellule saine, pourtant ce n'est pas du tout le cas.
La seule activité plus importante pour le modèle cancéreuse est celle de la phosphofructokinase. Or il s'agit d'une enzyme censée être sous-exprimée dans une cellule cancéreuse.
barplots_dict_2["pentose phosphate pathway"][0]
barplots_dict["pentose phosphate pathway"][0]
Le métabolisme du folate est supposé être suractivé dans une cellule cancéreuse.
C'est bien ce qu'on observe ici avec une activité beaucoup plus importante de DHFR, qui est d'ailleurs une cible thérapeutique anticancéreuse connue.
barplots_dict_2["folate metabolism"][0]
barplots_dict["folate metabolism"][0]
On s'attend à une activité très réduite des protéines de la beta-oxydation (ici, toutes les beta oxydations sont regroupées : mitochondriale et peroxysomale)
On observe effectivement une activité inexistante (flux nuls)
barplots_dict_2["beta oxidation"][0]
barplots_dict["beta oxidation"][0]
On s'attend à une activité très réduite des protéines associées à la navette carnitine.
On observe effectivement une activité très réduite de ces protéines.
barplots_dict_2["carnitine shuttle"][0]
barplots_dict["carnitine shuttle"][0]
Subsystem_fluxes_cancer = vu.build_reaction_df(HepG2, False)
treemap_dataframe = pd.DataFrame(Subsystem_fluxes_cancer["Glycolysis / Gluconeogenesis"])
title = "Glycolysis reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(treemap_dataframe, HepG2, title, path = ["compartment", "id"], color_by="compartment")
Subsystem_fluxes_healthy = vu.build_reaction_df(iHep, False)
treemap_dataframe = pd.DataFrame(Subsystem_fluxes_healthy["Glycolysis / Gluconeogenesis"])
title = "Glycolysis reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(treemap_dataframe, iHep, title, path = ["compartment", "id"], color_by="compartment")
title = "Cytoplasm reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_C,HepG2,title).show(renderer = 'notebook')
title = "Cytoplasm reactions -- <b>HEALTHY MODEL <br />"
vu.plot_treemap(df_Ci, iHep, title).show(renderer = "notebook")
barplots_dict["C_c"][0]
title="Mitochondrial reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(df_Mi, iHep, title).show(renderer = "notebook")
title="Mitochondrial reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_M, HepG2, title).show(renderer = "notebook")
barplots_dict["C_m"][0]
title="Peroxysomal reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(df_Pi, iHep, title).show(renderer = "notebook")
title="Peroxysomal reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_P, HepG2, title).show(renderer = "notebook")
barplots_dict["C_p"][0]